function [V,T] = mesh(nelx,nely,h)
% Generate an analysis domain and mesh. Nodes and Elements numbers are set as follows.
%
%11-12-13-14-15 %Node number
% | 5| 6| 7| 8| %Element number
% 6- 7- 8- 9-10 %Node number
% | 1| 2| 3| 4| %Element number
% 1- 2- 3- 4- 5 %Node number
%
V = zeros((nelx+1)*(nely+1),2); %Coordinates of each node
T = zeros(nelx*nely,4); %Index of nodes constructing each element
for j = 1:nely+1
    for i = 1:nelx+1
        V((j-1)*(nelx+1)+i,:) = [(i-1)*h (j-1)*h];
    end
end

% noisy the mesh
% for j = 2:nely
%     for i = 2:nelx        
%   %       V((j-1)*(nelx+1)+i,1) = V((j-1)*(nelx+1)+i,1) + rand*h*0.2;
%   %       V((j-1)*(nelx+1)+i,2) = V((j-1)*(nelx+1)+i,2) + rand*h*0.2;
%         V((j-1)*(nelx+1)+i,1) = 1 + sin(0.5*pi*(V((j-1)*(nelx+1)+i,1)-1));
%         V((j-1)*(nelx+1)+i,2) = (1 + sin(pi*(V((j-1)*(nelx+1)+i,2)-0.5)))/2;
%      end
% end

% originally, the node is not stored in an anti-clockwise terms!!
for j = 1:nely
    for i = 1:nelx
        T((j-1)*nelx+i,:) = [(j-1)*(nelx+1)+i (j-1)*(nelx+1)+i+1 j*(nelx+1)+i+1 j*(nelx+1)+i];
    end
end